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Breakdown of self -similar evolution in homogeneous perfect fluid collapse 
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The stability analysis of self-similar solutions is an important approach to confirm whether they act 
as an attractor in general non-self-similar gravitational collapse. Assuming that the collapsing matter 
is a perfect fluid with the equation of state P = up, we study spherically symmetric non-self-similar 
perturbations in homogeneous self-similar collapse described by the flat Friedmann solution. In the 
low pressure approximation k < 1, we analytically derive an infinite set of the normal modes and 
their growth (or decay) rate. The existence of one unstable normal mode is found to conclude that 
the self-similar behavior in homogeneous collapse of a sufficiently low pressure perfect fluid must 
terminate and a certain inhomogeneous density profile can develop with the lapse of time. 
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> ■ I. INTRODUCTION 

o ; 

Spherically symmetric self -similar gravitational collapse of a perfect fluid with pressure P given by the equation of 
state P — ctp is one of the most extensively-studied phenomena in general relativity. Many efforts have been made to 
solve the Einstein's equations governing its dynamics, which are reduced to a set of ordinary differential equations 
with respect to the single variable z = r/t. The flat Friedmann solution is well known as the unique analytically- 
found exact solution regular at the center and has played an important role in finding a family of solutions regular at 
the center [1, 2]. The homogeneous collapse described by this solution has been considered as the most basic process 
to spacelike singularity formation in the self-similar dynamics, while essential features of inhomogeneous collapse 
have been understood mainly through the detailed analysis of the general relativistic Larson-Penston solution [3] . 
In addition, it is also noteworthy that the perfect fluid critical collapse corresponding to the threshold of black hole 

■ formation has been confirmed to be described by one of the self-similar solutions [4, 5, 6]. 

\Q \ Several works have been devoted to the stability analysis of such self-similar solutions for spherically symmetric 
non-self-similar perturbations and have given important implications to more general non-self-similar gravitational 
""q I collapse. In particular, it is remarkable that the flat Friedmann solution and the general relativistic Larson-Penston 
Q" 1 solution were numerically confirmed to be able to act as an attractor in general spherically symmetric gravitational 
collapse for a. lying in the range < a < 0.036 [7]. In addition, the critical phenomena are illustrated in terms of the 
(3J[), time evolution of the single unstable normal mode which was found by the numerical analysis of the perturbations 
in the critical collapse [8, 9]. (See [10] for a recent review on the role of self-similar solutions as an attractor and the 
. i-h ', critical phenomena.) Although it may be interesting to study the stability problem more extensively, these numerical 

■ results should be confirmed through an analytical treatment of the perturbations. 

Recently, we have developed an analytical scheme to treat the stability problem by constructing the single wave 
equation governing non-self-similar spherically symmetric perturbations [11], which is reduced to the ordinary dif- 
ferential equation if we assume the perturbations to have the time dependence given by exp (ico log 1 1 \ ) . In this paper, 
using this analytical scheme, we study the stability problem for the flat Friedmann solution in the low pressure limit, 
i.e., < 0i <C 1. Fortunately, in the expansion with respect to the small parameter a, we can explicitly solve the 
master ordinary differential equation for the normal modes and consequently find the single unstable normal mode, 
which was not found in the numerical analysis [7]. 

We begin with a brief description of the perturbation theory for the flat Friedmann solution in Sec. II. In Sec. Ill, 
a discrete set of the normal modes and their growth (or decay) rate are derived in the low pressure limit, and the 
self-similar behavior turns out to be unstable in homogeneous collapse of a sufficiently low pressure perfect fluid. In 
Sec. IV, we see density inhomogeneities generated by the normal modes and explain the affection of the background 
transonic flow upon the growth (or decay) of the normal modes. In the final section, we summarize this paper and 
give a suggestion for the result of the numerical study [7]. In addition, we discuss an implication of the breakdown 
of the self-similar evolution in relation to critical phenomena. 
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II. PERTURBATION THEORY FOR SELF-SIMILAR HOMOGENEOUS PERFECT FLUID COLLAPSE 



In this section, we briefly illustrate our analytical scheme to treat spherically symmetric non-self-similar pertur- 
bations in homogeneous self-similar perfect fluid collapse described by the flat Friedmann solution (see [11] for the 
details). The line element considered throughout this paper is given by 

ds 2 = - e Mt,r) dt 2 + e 2MM) dr 2 + R 2 (f/ r) (^q2 + gin 2 (21) 

with the comoving coordinates t and r. In addition, the energy momentum tensor of a perfect fluid is expressed as 

T ab = (p + P)u a u b + Pg ab , (2.2) 

where p, P and the vector u" are energy density, pressure and fluid four velocity, respectively. As was mentioned in 
Sec. I, its equation of state is assumed to be 

P = otp (2.3) 

with a constant a lying in the range < a < 1. To discuss the self-similar behavior later, we use a new variable z 
defined by z = r It, instead of r. In addition, instead of the four unknown functions v, A, R and p, we introduce the 
following dimensionless functions: 

S[t,r) = -, rj(t,r) = 8Tir 2 p , M(t,r) = — , V(t,r) = ze x ~ v , (2.4) 

where the function m(t,r) is the Misner-Sharp mass. The function V is interpreted as the velocity of a z — const 
surface relative to the fluid element. 

From the Einstein's field equations for the system (2.1), (2.2) and (2.3), we can obtain the four partial differential 
equations governing the functions v, A, S and rj. By the two equations, the metrics v and A are explicitly given by the 
functions S and rj. The relations between the functions S and rj given by the remaining two equations become the 
simpler first order partial differential equations if we use the function M, which is explicitly given by the functions 
S and rj. Therefore we hereafter focus our concern on the functions S, rj and M. 

Now we consider spherically symmetric non-self-similar perturbations in the flat Friedmann background by ex- 
pressing the solutions for the Einstein's equations as 

S(f,z) = S B (z) {l + eSi(£,z)+0(e 2 )} , r]{t,z) = rj B (z) {l + erj X (i,z) + 0(e 2 )} , 

M(f,z) = M„(z) {l + eM 1 (t,z) + 0(e 2 )} (2.5) 

with a small parameter e. The functions Sg, rjB an d Mb for the flat Friedmann solution are given by 

S B (z)oc(-z)-P, // B (z)«z 2 , M b (z) « {-zf-^V (2.6) 



with the constant p defined as 



Note that the line element (2.1) for the flat Friedmann solution can be reduced to a form more familiar in cosmology 
through the coordinate transformation r^?exr 1- P[l]. 

As seen from the behaviors such that Sb — > and tjg —> oo in the limit z —> oo, a big crunch singularity appears at 
t = in the flat Friedmann background. Therefore we hereafter consider time evolution of the perturbations during 
t < (i.e., z < 0). As was mentioned in Sec. I, in such a region there are two characteristic surfaces at which a 
constraint is imposed on the behavior of the perturbations. One is the regular center located at z — 0~, and the other 
is the sonic point located at z = z s defined as a point at which the velocity of a z = const surface relative to the fluid 
is equal to the sound speed, i.e., V B = — ^foc, where the function Vg for the background flat Friedmann solution is 
found to be 



y B ( z )oc(-z) 1 -p 



(2.8) 
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(2.13) 



Using the perturbation equations for Si, t]i and Mi, we can obtain the following single wave equation: 

Xuu ~ % + W(0%u + F(£)Y£ + = (2.9) 

for the function Y defined as 

Y(f,z) = S!(f,z) - /(z)Mi(f,z) , (2.10) 

where the functions W, F, U and / are given in Appendix A. The new variables u and £ are related to the variables f 
and z by 

« = log(-f) + I(z), i( 2 ) = 5f J-^log(l- )t 2 (z)), (2.11) 
where the function x is defined as 

X = . 

The remaining perturbation equations yield the relations 

Mi(f,z) = Bi(z)Y(t,z) + B 2 (z)Y'(f 7 z) + B 3 (z)Y(f / z) / (2.14) 

7i(f 7 z) = M x (t,z) + — 1 , M' 1 (t,z)-3S 1 (t,z) - -}—S[{t,z) , (2.15) 
3(1 -p) 1-p 

where the dot and the prime represent the partial derivative with respect to log |f | and log |z|, respectively, and 
the functions B\, B 2 and B 3 are also given in Appendix A. The perturbations Si, r\\ and Mi are determined by the 
solution Y for the wave equation (2.9) via Eqs. (2.14), (2.10) and (2.15). 
The wave equation (2.9) allows us to consider the modes (p defined as 

Y(f,z) = <p(x,co)e iw{ - u+ ® (2.16) 

with the spectral parameter co. It is mathematically convenient to use the variable x, instead of £, which can cover 
the whole region between the regular center and the sonic point in the finite range < x < 1, irrespective of the 
parameter value a. Then the equation for <p is found to be 

2ico — 2(1 — p)x — F ico(F + W) + U „ 

*** + (i-p)d -#) *> ~ (i-mi-*)> * = • (217) 

We will obtain the normal modes as the solutions for Eq. (2.17) satisfying the boundary conditions such that cp is 
analytic both at the regular center x = and at the sonic point x — 1. (The same boundary conditions were required 
in the numerical analysis done by [7].) The leading behavior of <p near x — 1 is given by 

</>~Ci(a;) + C 2 (aO(l-x) fc , (2.18) 

where k — {(l + z'a;)/(l — p) } + 1 . The ratio C 2 / C\ will be uniquely determined by the requirement of the analy ticity 
of (p at x = 0, and discrete eigenvalues to — lo„ giving the normal modes <p(x,Lo n ) = <p n (x) will be derived by the 
equation 

C 2 (o;) = 0. (2.19) 
Here we would like to note that there exists an exact solution <p = (p s for Eq. (2.17) written by 

r fi _i_ r \3(l-«)/(3«+l) 
rg 2x + l+3a v ' 

if the spectral parameter co is equal to co g defined as 

Wo = . (2.21) 

s 1 + a 

This solution <p s is one of the normal modes (p„ but corresponds to a gauge mode due to an infinitesimal transforma- 
tion of t. In fact, all the perturbations r\\, S\ and Mi obtained from <p g are found to be independent of z. Although 
the gauge mode is obviously unphysical, the presence of such an exact solution will be mathematically useful for 
checking the validity of the analysis of Eq. (2.17). 
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III. NORMAL MODES IN THE LOW PRESSURE LIMIT 



Although the master equation (2.17) for the perturbations in the flat Friedmann background is given by a simpler 
form (if compared with the form for any other self-similar backgrounds), it is still a difficult task to solve the eigen- 
value problem, and in [11] the absence of the unstable normal modes was clearly proven only in a limited range of 
to. Such a difficulty may be overcome if we consider the low pressure limit a. — > 0, keeping the variable x finite in 
the range < x < 1 and expanding the solution <p(x, to, a) analytic at x = as follows, 

CO 

(p(x, to, a.) = Yj <* i( P® {x,to) . (3.1) 
i=0 

For the lowest-order solution cp^ we obtain the equation 

Ccp^\x,to) =0 , (3.2) 
where the ordinary differential operator £ is given by 

d 2 3 f . 2x 3 +3x 2 + 10x + 3l d 
C = + - ~ ^ 2ito + 



(3.3) 



dx 2 1 -x 2 \ 3x(l + 2x) J dx 

2x 4 + 4x 3 - x 2 - 2x - 1 



' x(l+2x)(l + x)(l-x 2 ) \ 



2ito(2x 2 -2x- 1) 



We would like to emphasize that this limit does not mean to consider an exactly pressureless fluid (i.e., a dust fluid) 
because the requirement of the analyticity of (p at the sonic point x — 1 is not missed. The crucial point in this 
approach is that we can explicitly derive general solutions for Eq. (3.2). In particular, the solution <p°) satisfying the 
boundary condition at x — is written as 

(f>( \x,to) = Zt(x,to) - Z 2 {x,to) , (3.4) 

where the functions Z\ and Z2 are the two independent solutions for Eq. (3.2) and given by 

(l+x) 3 {-6x 2 to 2 -4x(x 2 -3x + l)ito + (l-x) i } , 

Zi(x,co) = ± i-k 3' ^ r U -, 3.5 

x^yl + 2x) 

^ . s (l-x) i+3 ^{l + x) 1 - 3lw {x 2 + 2{l + icv)x + l} 

Z2{X ' W) = x\l + 2x) ' (3 - 6) 

It is clear that this solution ^>(°) becomes analytic also at the sonic point x = 1 if the spectral parameter to is equal 
to co^ given by 

4 0) = (3.7) 

with n defined as non-negative integers. However, for n — 0, 2 and 4, the function <p^ turns out to vanish. Hence, 
the values of a^ ' are given only for n = 1, 3, 5, 6, ■ ■ ■ . It can be easily found that the value of tof 1 ^ and the function 
are identical with the value of to g and the gauge mode <p g in the limit a — > 0. While this assures that 

to^ represents the eigenvalues approximately written in the limit a — > 0, it may be unclear how the condition (2.19) 
to obtain the eigenvalues to n is used in this approach. 

To discuss this point, we consider the next order solution <p( l \x,to) in Eq. (3.1), using Eq. (3.4). The inhomoge- 
neous ordinary differential equation for the function <p^> is written as 

Ccp^(x,to) = ](x,to) , (3.8) 

where the function / is given in Appendix A. Here we note that the imaginary part of co^ is smaller than 4/3. 
Because we consider the solution analytic at x = for to nearly equal to to^\ from Eq. (3.8) we obtain 

*V>{xM = Wxrt f ^ ){y '^' w) dy + ^{x,cv) ^ z 2 (y,c,)7(y^) (0) 

Jo w(y,to) Jx w(y,to) 
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where a(co) is an arbitrary constant and the Wronskian w of (p^ and Z 2 is given by 

8(1 - x)^ l+ico \\ + x)^ l ~ iw ^co(co - 2i) (3a; - 2z) (3a; - 4i) 



(3.10) 



x 3 (l+2x) 2 

The eigenvalues cv n (a) giving the normal modes <p„ up to the first order of a will have the form 

Lo n {d) = o4 0) +o4 1) a + 0(a 2 ) . (3.11) 
Because of this expansion of co„ (a), Eq. (3.1) for the normal modes <pn{= <p{x, co„ (a), a) ) can be rewritten to the form 



<Pn(x f *) = <$\x) + I ^(x,cor)^n'+<p K n'(x) J a + 0(a z ) , 



(3.12) 



where <pn\x) = tp^ l \x,co^). In this expansion scheme, the function <p n may be non-analytic at x = 1 owing to 
terms containing the logarithmic factor log(l — x) in 

8^(°)/aa;(x / 4 0) ) and It is straightforward to derive such 

non-analytic terms, and we have 



^-^(x,co^) — K„(x) - 3iZ 2 (x,co^) log(l - x) , 
<p)t\x) = L n (x) - b n Z 2 (x,co^)log(l - x) , 



(3.13) 
(3.14) 



where the functions K n and L n are analytic at x = 1 and the coefficient b n is obtained from the term proportional 

to (1 — x)~ x in the expansion of the function Z\] /w for co = co^ around x — 1. (Note that the function Z-yJ/w 
appears in the integrand of the first integral in the right hand side of Eq. (3.9) because of Eq. (3.4).) It is clear that the 

non-analyticity of <p n at x = 1 can be removed if co^ is chosen as follows, 

(1) _ b n i 



(3.15) 



To estimate the coefficient b n , we rewrite the function Z\J/w for co = co^ into the form 

B(x,4 0) ) {(1 - *) 3 -"M*,4 0) ) + (1 - x)-"/ 2 (x,4 0) ) + (1 - x)- l h{x,co^)} , (3.16) 



Zi(x,o;Ji 0) )/(x,a;i 0) ) 



1 (o)^ 
w{x,co„ ') 

where the functions /1, 72 and 73 are given in Appendix A and the function B is defined as 

Z 1 (x,co) 



B{x,co) 



w{x,o}){\- x)-^ +io} ) ' 
,(0) 



(3.17) 



Note that the functions B, 7i, } 2 and 73 for co = co n ' become finite at x — 1. From the form given by Eq. (3.16), we 
can easily find 



b n = B{l,co$ ) )h{l,uT)+bn, 



(OK 



where b n is given by 



h = B(l,o;f )7 2 (l,o;i 0) ), 



£3 = \ (B(x,u,W)J 2 (x,a,i 0) j) x 



x=l 



and for n > 5 

(_l)n-4 d n-4 



(n-4)! dx"" 4 



B(x,4 0) )7i(^4 0) )) ^ +^^f ^(B^^Mx,^)) 



n— 1 jn— 1 



(3.18) 

(3.19) 
(3.20) 

(3.21) 
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We can easily calculate the first term in the right hand side of Eq. (3.18) to be 

B(l,4 0) )/3(1,4 0) ) = -2(n - 1) . (3.22) 
In addition we find that the value of b„ is given by the unified form 

b» = 6(-l)« (3.23) 

for n — 1, 3, 5, 6, • • • . Thus we arrive at the result 

«»(«) = + - 1) + 2(-l) M J oci + 0(a 2 ) (3.24) 

for the eigenvalue problem Eq. (2.19) giving the normal modes <p„ written as 



<p„{x,ct) = (pli\x) + 



Kn(x) |-|(n - 1) +2(-l)"j i + L„(x) 



a + 0(a 2 ), (3.25) 



where the functions K n and L« are explicitly given through Eqs. (3.13) and (3.14). It should be noted that the value 
of lo\ given by the above equation is identical with the value of co g up to the first order of a, as was expected. 
Although the normal modes <p n appear to have an ambiguity due to the existence of the constant a in Eq. (3.9), it will 
be uniquely determined if we require the analyticity at x — 1 up to the second order of a. In fact, we can confirm 
that Eq. (3.25) for n = 1 is identical with the gauge mode (p g only for a = 7 — 12 log 2. 

The most important result is that the imaginary part of o;^ ' given by Eq. (3.7) is positive, namely, there exists one 
unstable normal mode, at least for sufficiently small values of a. The proof concerning the absence of unstable normal 

modes shown in [11] is not applicable to the normal mode cps with the small growth rate \m(w^ + aco^ ) obtained 
here. Hence the flat Friedmann solution does not act as an attractor in homogeneous collapse of a sufficiently low 

(1) 

pressure perfect fluid. However, the first order correction co^ has a negative imaginary part. This means that the 
growth rate of the unstable normal mode becomes smaller as the value of a. increases, and the self -similar behavior 
might be stable in homogeneous collapse of a higher pressure perfect fluid such as the radiation fluid (a = 1/3) 
and the stiff fluid (a = 1). This is an interesting problem to be further studied in future works. Moreover, it is also 
interesting that there exists an infinite set of the stable normal modes (i.e., (p„ for n > 5). In the next section, we will 
derive the density perturbation rji corresponding to the normal modes <p„ to see what configuration of perturbed 
inhomogeneous fields can develop or decay with the lapse of time. 

IV. DENSITY PERTURBATIONS 

Let us denote the perturbed density r\\ corresponding to the normal modes cp„ by r]\( n y which is expanded with 
respect to a as follows, 

7 1(n) (f, x, «) = / 7 ^ ) (x)exp{/4 0) log(-0}+O(a) . (4.1) 

Here we focus our attention on the leading term depending on x. Through Eqs. (3.4), (2.16), (2.11), (2.12), (2.14), 
(2.10) and (2.15) in the limit a — > with the eigenvalues w = u>n\ we obtain 

= ^[(l-*)"~ 2 {3* 2 + 3(«-2)*+(n-l)( n -3)} 

-(l+x)"- 2 {3x 2 -3(n-2)x+(n-l)(«-3)}] . (4.2) 

In Fig. 1, we show the configuration of the density perturbation j/j?^ (x) normalized by its value at x = 0, namely, 

<>> = - ( "~ 4 ""~3 2) " ( " +2> - <«) 
for n = 1,3,5,6,7,8,9 and 10 in the range < x < 1 . It is shown in this figure that the normal mode 
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corresponds to the gauge mode r]g = const at a given time as was mentioned in Sec. II. In addition, it should be 
emphasized that the normal mode rj-f^y which is also constant at any x, is a physical normal mode because the 
corresponding perturbation Mj depends on x in the range < x < 1. 

From Fig. 1, we note that the amplitude of all the stable normal modes (i.e., the normal modes for n > 5) at the 
sonic point x = 1 remains non-zero and the amplitude of the normal modes for n > 7 rather increases towards the 

sonic point x = 1 from the center x = 0. The ratio vf^) (1) /^(n) ^ mcreases as n becomes larger and the decay rate 
— Im(o; n ) given by Eq. (3.24) has the same tendency. This implies that the density perturbation generated near the 
sonic point is rapidly carried away to the supersonic region x > 1 by the background transonic flow and the growth 
of such a density perturbation in the subsonic region is prevented. It is remarkable that only for n = 3, the value of 

7^ vanishes at x = 1. This seems to be the most favorable configuration of rji to allow the growth of the density 
perturbation due to the effect of its own self -gravitation. Further, Eq. (3.24) clearly shows that the imaginary part of 
co„ decreases as a. increases. The dispersive effect due to the pressure against the self-gravitation can enhance the 
above-mentioned decay process of the perturbations. 

In this section, we have focused our concern on the configuration of the normal modes in the subsonic region 
< x < 1. This is mainly because the eigenvalue problem was set under the boundary conditions at x = and 
x = 1 and from the viewpoint of causality, any disturbances in the supersonic region cannot affect the process in 
the subsonic region. Further, as another reason, we would like to point out that the approximation a ~ to derive 
Eq. (3.24) becomes mathematically unreliable in the region far away from the sonic point, i.e., x ^> 1. If one tries to 
understand more global features of the normal modes by analyzing Eq. (2.17) in the limit a — > 0, the terms included 
in Eq. (2.17) which can be negligible in the subsonic region but become important for x 2> 1 should be taken into 
consideration. 



V. SUMMARY AND DISCUSSION 



In this paper, we have studied the stability problem for the self-similar behavior in homogeneous collapse of a 
perfect fluid with the equation of state P — ocp, using the perturbation theory developed in [11]. We have derived 
the single ordinary differential equation (2.17) governing spherically symmetric non-self-similar perturbations with 
the time dependence exp {/a? log ( — t)} in the flat Friedmann background and set up the eigenvalue problem to 
determine the value of the spectral parameter co. In the low pressure approximation a — > 0, we have succeeded in 
deriving explicitly an infinite set of the eigenvalues and the normal modes given by Eqs. (3.24) and (3.25). Because 
one of such normal modes is an unstable normal mode, we have concluded that non-self-similar inhomogeneous 
disturbances can develop in homogeneous collapse of a sufficiently low pressure perfect fluid. 

As was mentioned in Sec. I, the unstable normal mode obtained in this paper was not found in the numerical 
analysis [7], which has rather claimed that for a lying in the range < a < 0.036, the geometrical structure and the 
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fluid motion at late stages in general non-self -similar collapse of a perfect fluid (which is initially homogeneous) can 
be well described by the flat Friedmann solution. However, recalling that the unstable normal mode obtained in this 
paper has the small growth rate (i.e., Im(o?3)) less than 1/3, we suggest that even if the geometrical structure and 
the fluid motion once become similar to those of the flat Friedmann solution, non-self -similar disturbances become 
significant at much later stages which were missed in the numerical simulation [7]. Although the general relativistic 
version of the Larson-Penston solution was also suggested to act as an attractor in general inhomogeneous collapse 
by the results of the numerical simulation and the normal mode analysis in [7], we claim that its stability should be 
also confirmed in our analytical scheme. 

Because of the above-mentioned result of the numerical simulation [7], the result obtained in this paper allows us 
to interpret the flat Friedmann solution as an intermediate attractor in general non-self -similar perfect fluid collapse 
which starts from a nearly homogeneous density profile, at least for sufficiently small a. If the general relativistic 
Larson-Penston solution is confirmed to be stable, a transition from the flat Friedmann stage to the general relativistic 
Larson-Penston stage may occur in gravitational collapse. (Such a transition was also mentioned in [7].) The so-called 
critical solution corresponding to the threshold between the black hole formation and the complete dispersion of the 
fluid is the well-studied self-similar perfect fluid solution (see e.g., [6]) acting as an intermediate attractor in general 
inhomogeneous collapse. It is interesting to note that there is a common feature between the flat Friedmann solution 
and the critical solution, namely, the single unstable normal mode exists for these solutions. It was found from the 
idea of the renormalization group in [8, 9] that such a feature is essential to the scaling-law and the universality 
observed in the critical phenomena. Therefore what critical phenomena are relevant to the flat Friedmann solution 
will be an interesting problem to be investigated in future works. 



APPENDIX A: FUNCTIONS IN THE PERTURBATION EQUATIONS 



It is an easy task to calculate the functions involved in Eqs. (2.9) and (2.14) by using the formulae given in [11]. 
The results are summarized as follows, 

(1 + 3a)x 3 + (7 + 9a)x 2 + (1 + 3a) (7 + 9a)x + (1 + 3a) 2 

3(l + a)x(l + 3ft + x) ' ( ' 

6(1 - a)x 4 - (1 + 3ft) (9ft - ll)x 3 + 3(1 + 3ft)(5 + 7a)x 2 + (1 + 3a) 2 (13 + 3a)x + 3(1 + 3a 3 ) 

/ (A2) 



U 



3(1 + a)x(l + 3ft + x) (1 + 3a + 2x) 
(1 + 3ft) (1 - x 2 ) {2(ft - l)x 4 + 4(ft - 1) (1 + 3a)x 3 - (1 + 3a) (5a - l)x 2 + 2(1 + 3a) 2 x + (1 + 3a) 3 } 

3(1 + a) 2 x 2 (l + 3a + x)(l + 3a + 2xj 



_ 2ax - 3a - 1 

1 3ft(l + 3ft + 2x)' K ] 

9a(l + 3a + 2x) 2 

1 8(l + 3ft)(l + 3ft + x) ' K ] 
= 9a(l+x)(l+3a + 2x) 2 

2 8(l+3a)x(l + 3a + x) ' 1 ' 
9ft(l + 3a + 2x) {2(a - l)x 2 + (1 + 3a) 2 } 

3 ~~ 8(l + a)(l+3a)x(l+3a + x) ' ( ) 

If the functions involved in Eq. (2.17) are expanded with respect to ft as follows, 

F + 2(1 - p)x = Q (x) + Qi (x)a + 0(a 2 ) , (A7) 
F + W = P (x)+P 1 {x)a + O(oc 2 ) , (A8) 
U=U {x) + U 1 (x)a + 0{oc 2 ) , (A9) 

the inhomogeneous term / in Eq. (3.8) is given by 



/(-v,u'l - ~y^2 {-6(2ia; - Q ) - 3Q 1 } ^ + _* {-36(fa;F^ + U )+ 9{icoP r + Uj)} ^ , (A10) 
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which leads us to the final result 

}(x,co) 



(l-x)(l + x) 2 +3^X 4 (l+2x) 4 

x (1 + x) 3+3iw | (x - l) 3 (64x 6 + 122x 5 + 119x 4 + 74x 3 + 34x 2 + 16x + 3) 

-2(104x 8 - 102x 7 - 45x 6 + Ax 5 - 109x 4 - 43x 3 + 2x 2 + 7x + 2)iw 

-2x(160x 6 - I22x 5 - 217x 4 + 67x 3 + 98x 2 + 28x + 4)cv 2 + 36x 4 (l + 2x) (3 + 4x)io; 3 } 

+ (1 - x) 3+3iw {(1 + x) 3 (64x 6 + 122x 5 + 119x 4 + 74x 3 + 34x 2 + 16x + 3) 
+2(88x 8 + 420x 7 + 880x 6 + 1047x 5 + 821x 4 + 447x 3 + 155x 2 + 28x + 2)iw 

-2x{l + 2x) (56x 5 + 238x 4 + 336x 3 + 225x 2 + 73x + 8)o; 2 - 24x 2 (l + x) (1 + 2x) 3 ^ 3 } . (All) 
In addition, the functions ]\, ]i and ]^ in Eq. (3.16) are given by 

h(x,u>) = -^^^4(64x 6 + 122x 5 + 119x 4 + 74x 3 +34x 2 + 16x + 3) / (A12) 

h(*,w) = t \\ I ~2(104x 8 - 102x 7 - 45x 6 + 4x 5 - 109x 4 - 43x 3 + 2x 2 + 7x + 2) 

x 4 (l + 2x) 4 L 

+2x(160x 6 - 122x 5 - 217x 4 + 67x 3 + 98x 2 + 28x + 4)iu> + 36x 4 (l + 2x) (3 + 4x)w 2 } , (A13) 



h(x,cv) 
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(l + x) 2 + 3i ^x 4 (l+2x) 4 

x | (1 + x) 3 (64x 6 + 122x 5 + 119x 4 + 74x 3 + 34x 2 + 16x + 3) 

+2(88x 8 + 420x 7 + 880x 6 + 1047x 5 + 821x 4 + 447x 3 + 155x 2 + 28x + 2)iw 

-2x(l + 2x)(56x 5 + 238x 4 + 336x 3 + 225x 2 + 73x + 8)o; 2 - 24x 2 (l +*)(! + 2x) 3 fo; 3 } . (A14) 
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